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PROGRAM SUMMARY 



Title of program: WWGENPV 

Computer: DEC VAX, HP/APOLLO, IBM/RS6000; Installation: INFN, 
Sezione di Pavia, via A. Bassi 6, 27100 Pavia, Italy 

Operating system: VMS, UNIX 

Programming language used: FORTRAN 77 

Memory required to execute with typical data: 500 KByte 

No. of bits in a word: 32 

The code has not been vectorized 

Subprograms used: random number generator RANLUX [1]; routines from the 
CERN Program Library are also used. 

No. of lines in distributed program, including test data, etc.: 2135. 

Correspondence to: 
MONTAGNA@PV.INFN.IT; 

NICROSINI@VXCERN.CERN.CH, NICROSINI@PVINFN.IT; 
PICCININI@PV.INFN.IT 

Keywords: high energy electron-positron collisions, VF-pair off-shell produc- 
tion, four-fermion final state, LEP2, QED corrections, electron structure func- 
tions, experimental cuts, Monte Carlo integration, event generation, impor- 
tance sampling. 

Nature of physical problem 

The precise measurement of the FF-boson mass M w constitutes a primary 
task of the forthcoming experiments at the high energy electron-positron col- 
lider LEP2 (2Mw < — 210 GeV). This experimental goal will be possible 
through the study of the four-fermion reaction e + e~ — > W + W~ — > 4f. A 
meaningful comparison between theory and experiment requires an accurate 
description of the fully exclusive process e + e~ — > W + W~ — > 4/, including the 
main effects of radiative corrections, with the final goal of providing predictions 
for the realistic distributions measured by the experiments [2]. 

Method of solution 

In order to efficiently perform the high-dimensional integration (due to the 
four-fermion final state) in the presence of realistic cuts, the total cross section, 
the energy and invariant-mass loss and the W invariant-mass distributions are 



1 



computed by means of a Monte Carlo integration for weighted events. For 
simulation purposes, the program can also be used as an event generator that 
provides a sample of unweighted events, defined as the components of the 
four final-state fermions momenta, plus the radiative variables of the incoming 
electron and positron, plus yfs, stored into proper n-tuples. The importance- 
sampling technique [3] is employed to take care of the peaking behaviour of 
the integrand, both in the integration and in the generation modes. 

Restrictions on the complexity of the problem 

Electroweak processes originating the same four-fermion final state of the 
double-resonating diagrams e + e _ — > W + W~ — > 4/ {background processes) 
are not included in the present version of the program. The treatment of the 
initial-state radiation is limited to the gauge-invariant subset of the leading- 
logarithmic corrections which are the most important contribution from the 
phenomenological point of view. Furthermore, the kinematics of the event in 
the presence of the initial-state radiation is limited to the collinear approxi- 
mation, i.e. Pt/pl effects are neglected for the time being. 

Typical running time 

As integrator, the code needs about 5 h of HP 9000/735 for generating 10 8 
weighted events; the relative errors on the cross section and on the energy 
and invariant-mass losses are about 2 x 10~ 4 and 6 x 10~ 4 , respectively. The 
generation of a sample of 10 4 unweighted events requires about 3 min on the 
same system. 

Unusual features of the program: none 
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LONG WRITE - UP 



1 Introduction 

The precise measurement of the W^-boson mass M^, and more generally the 
study of its physical properties, constitutes a primary task of the forthcoming 
experiments at the high energy electron-positron collider LEP2 (2M\y < ^/s < 
210 GeV). 

A large amount of theoretical work has been performed up to now, following 
essentially two main streams. 

On the one hand, one can find semi-analytical works, whose leitmotiv is 
computing analytically as much phase-space integrations as possible |IJ. Typ- 
ically, these works give (QED-corrected) distributions for the invariant masses 
of the W^, integrating a priori all the other independent phase-space vari- 
ables, from which total cross sections can be easily computed. After such an 
analytical effort, the semi-analytical codes developed along this line are very 
fast, precise and efficient. They are hence typically well suited for fitting pur- 
poses, their main drawback being the lack of a full reconstruction of the true 
event phase-space configuration. 

On the other hand, one can find Monte Carlo codes. In principle they can 
guarantee the full event reconstruction, since they handle the four-momenta 
of the final-state particles. Of course, the accuracy with which this goal is 
achieved depends on the dynamics built-in in the code; in particular, the ex- 
clusive control of all the features of the full four-fermion event, including all 
the possible correlations between final-state particles, requires the knowledge 
of the exact four-fermion matrix element, at least in the Born approximation 
(for instance, no isotropic decays of the intermediate bosons can be allowed). 

Common to both the above-mentioned approaches, is the problem of treat- 
ing QED corrections, which, as is well known, are large and strongly dependent 
on the experimental cuts applied. In particular, it happens that the energy lost 
in the beam pipe (which at the leading-logarithmic level is well described by 
initial-state radiation in the collinear approximation) represents an important 
part of the systematic error on Mw- Therefore one needs also a reliable model 
for, at least, initial-state radiation. 

Some Monte Carlo codes are already available for the four-fermion pro- 
cesses of interest at LEP2 0. Nevertheless, the critical comparison between 
independent approaches to a problem has always been a useful and powerful 
tool for understanding the theoretical accuracy of the results. The aim of this 
work is to present and describe in some detail WWGENPV, a Monte Carlo event 
generator for e + e _ — > W + W~ — > 4/, including initial-state QED corrections 
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in the structure functions formalism. 



2 Theoretical Background and General Fea- 
tures 



The theoretical formulation implemented in the program has been described in 
some detail in ||. To calculate the cross section of the e + e~ — > W + W~ — > 4/ 
process in the presence of the initial-state radiation, the basic formula reads: 

a(s) = J dx x dx 2 D(x h s)D(x 2 , s)d[PS) jj^fL ■ (!) 

In this equation, d[PS] denotes the volume element in the 7-dimensional phase 
space connected to the independent variables of the hard scattering reaction in 
the centre-of-mass frame; dao/d[PS] is the exact differential kernel cross sec- 
tion in the Born approximation, limited at present to the P^-signal diagrams. 
The inclusion of the complete set of background diagrams is in progress. At 
any rate their effect can be cut-off by proper selection criteria on the final-state 
invariant masses. Finally D(x, Q 2 ) is the electron structure function describ- 
ing initial-state radiation ||]. For its explicit expression the following form is 
chosen: 

exp{|/3 (f-7fi)} P g_! P 

D(x, s) = ^-y — v i ^ — — — (1 — X) 2 — —(1 + x) 

In T 

(2) 



r (i + \p) 



+-I3 2 



-4(1 + x) ln(l - x) + 3(1 +x)lnx- 4-^- - 5 - x 

1 — x 



with 



P = 2-{L-\), (3) 

7T 

where L = In (s/m 2 ) is the collinear logarithm, 7^ is the Euler constant and 
T(z) is the gamma function. The first exponentiated term of Gribov-Lipatov 
form describes multiphoton soft radiation, the second and third ones hard 
collinear bremsstrahlung. This choice of the structure function is a consis- 
tent solution of the Renormalization Group evolution equation at the leading- 
logarithmic level, and guarantees that, after an expansion down to 0(a), the 
well-known factorized and hence universal bremsstrahlung spectrum is exactly 
recovered. 
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With the choice of the independent variables adopted in || , the differential 
kernel cross section is proportional to the tree-level squared matrix element as 
follows: 



dcr 
d[PS] 



oc 



(af4 + b/j,i + c)(a'r 2 2 + b'r 2 + d 



(4) 



Therefore, in the kernel cross section two kinds of "singularities" can be rec- 
ognized, namely the one related to the vector boson propagators (on fi 2 2 ), the 
second arising from the Jacobian appearing in (^) (on /i 2 ^). To efficiently 
perform the Monte Carlo integration in the presence of the above peaking struc- 
ture, the importance-sampling technique || has been employed. For instance, 
for the M^-boson propagator singularities, an integration must be performed, 
which has the form 



f(x 2 



,2^2 



+ b 2 ' 



(5) 



where x 2 = /if, a 2 = M^/s, b 2 = T 2 M^/s 2 and the function f(x 2 ) is regular. 
Multiplying and dividing the integrand by the weight function 



p(x z 



N 1 (x 2 n} x\ I ) 
(x 2 ~a 2 ) 2 + b 2 ' 



(6) 



with 



N{x 2 m ,x 2 M ) 



dx 



r 



,2\2 



one obtains 



Arctg 



dx 2 p(x 2 



x 



M 



+ b 2 

— Arctg 



(7) 



P l (x 2 )f(x 2 



X' 



,2\2 



+ b 2 



where the quantity in square brackets is no longer singular at x 2 
integral / can be rewritten as 

Jo (x 2 (r) - a 
introducing the cumulative variable c as 



p-\x 2 (c))f(x 2 (c)) 

2\2 _|_ Jj2 ' 



c{x 2 ) 



dz 2 p{z 2 ). 



The 



(9) 



(10) 
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Therefore, by generating c flat between and 1 an accurate estimate of / can 
be obtained. 

Concerning the singularities arising from the Jacobian factor, the integra- 
tions to be performed are of the following kind: 

," M dy , SM (11) 

Vm v (2/ - 2/1X2/2 - y) 

In this case the importance sampling can be applied by introducing the weight 
function 

p( y ) = N- 1 {y m ,y M ){y2 - y)~*(v - yi)~K ( 12 ) 

with 

rVM _i _i 

N(y m ,y M ) = / dy(y-y m ) *{y M -y) 2 
•>ym 



Arcsinw — — — — Arcsin< 
V 2/2 - Vi 





-l/i 


2/2 


-2/1 



(13) 



When including initial-state QED corrections, in order to perform the im- 
portance sampling on the generation of the radiative variable x in the soft 
region x — > 1, the weight function is chosen as the all-order soft part of the 
structure function, i.e. 

p{x) =N- 1 (x m ,x M )^(l-x)%-\ (14) 
with the normalization factor N(x m , %) given by: 

f x M fi p 

N(x m ,x M ) = / dy-(l-x)2 

= (l-x m )f -(l- XM )§. (15) 



With respect to the original formulation of QED corrections given in 0, 
besides a different choice of the structure function, the Coulomb singularity, 
which is important in the threshold region, is included in the code according 
to the very recent formulae given in 0. For details and links with other 
approaches to the problem, the reader is referred to the original literature. 



3 Program Structure 

After the initialization of the Standard Model parameters and of the electro- 
magnetic quantities, the calculation within the program proceeds as described 
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in the following. As a first step, the random number generator RANLUX is called 
to extract a vector of twelve components which are converted from REAL*4 to 
REAL*8 precision in order to be consistent with the calculational environment, 
which is in REAL*8 precision. The first two random numbers are used to gen- 
erate the radiative variables x\$. according to the distribution given by the 
soft exponentiated part of the structure function; eight random numbers are 
employed to fully reconstruct the kinematics of the four-fermion hard scatter- 
ing process in the centre-of-mass frame: seven of them are needed to generate 
the independent variables previously introduced for the kernel reaction, one to 
choose between the two possible solutions of the non-linear constraint. Given 
the independent variables, all the invariants related to the 2 — » 4 kernel re- 
action are computed. The variables //? )2 ,3 an d T2 are, in particular, generated 
by means of the importance-sampling procedure to take care of the peaking 
behaviour of the integrand, as described in the previous section. The other 
variables (£1,3,71) are generated according to the uniform distribution in the 
physical range relative to an extrapolated set-up. When initial-state radiation 
is switched on, the four-momenta of the outgoing fermions are reconstructed in 
the laboratory frame, taking into account the Lorentz boost due to the emis- 
sion of unbalanced electron and positron radiation and given by the (3l Lorentz 
factor (3l = {x\ — X2)/(xi + x%). At this point, one further random number is 
generated to perform a random rotation around the beam axis to account for 
a trivial overall azimuthal angle. The last random number is required in the 
generation branch for hit-or-miss. After that, a call to SUBROUTINE CUTUSER 
is performed in order to implement the rejection algorithm on the final-state 
particles. If the event is accepted, the squared matrix element in Born approx- 
imation is returned by a call to SUBROUTINE BORN. In conclusion, at the end of 
the loop over the events, the averages and squared averages are cumulated to 
provide the "weighted events" Monte Carlo results for the required observables 
and their numerical errors. 

When the program runs as a generator of unweighted events, after the hit- 
or-miss procedure, the components of the four final-state fermions momenta, 
plus the radiative variables of the incoming electron and positron, plus y/s, 
are stored into a proper n-tuple. Each event stored in the n-tuple has the 
structure (xi, X2, -E&, <fi, <?2, Q3, Q4), where Ej, is the beam energy and cfi are the 
three-momenta of the outgoing particles, assumed to be ultrarelativistic. In 
this case, the cross section for unweighted events is given as output. 
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4 Input 



Some input parameters and flags are introduced in the program in order to 
allow the user to select the output of interest. A short explanation of their 
meaning is addressed in the following and is also printed when running the 
code interactively. 

0GEN( CHARACTER* 1) 

It controls the use of the program as a Monte Carlo integrator for weighted 
events (OGEN = I) or as a Monte Carlo event generator of unweighted events 
(OGEN = G). If the first choice is requested, the code returns in a LOGFILE 
the value of the cross section for weighted events, together with the aver- 
age energy loss and invariant-mass loss, defined as ((2 — x\ — x 2 )\^s/2) and 
((1 — XiX2)\fs/2), respectively. If the second choice is requested, it prints 
in a LOGFILE the value of the cross section for unweighted events and fills 
an n-tuple, EXAMPLE.DAT, containing the generated (unweighted) events. The 
LOGFILE contains also the value of NBIAS/NMAX. If it is > 0, the generation 
procedure is biased and it is necessary to run the code again, after having 
properly increased the value of FMAX (default = 1000). 

RS(REAL*8) 

The centre-of-mass energy (in GeV). 
NHITWMAX (INTEGER) 

Required by the integration branch. It gives the maximum number of calls for 
the Monte Carlo loop. 

NHITMAX ( INTEGER) 

Required by the event-generation branch. It is the maximum number of hits 
for the hit-or-miss procedure. 

IQED (INTEGER) 

This flag allows the user to switch on/off the contribution of the initial-state 
radiation. If IQED = the distributions are computed in lowest-order approx- 
imation, while for IQED = 1 the initial-state QED corrections are included in 
the calculation. This choice can be performed both for the integration and for 
the generation mode. 
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ODIS (CHARACTER*!) 



Required by the integration branch. It selects the kind of experimental dis- 
tribution. For ODIS = T the program computes the total cross section (in 
pb) of the process, eventually including the cuts supplied by the user in 
SUBROUTINE CUTUSER (see below); for ODIS = W (only possible with OGEN = 
I) the W invariant-mass distribution is returned (in pb/GeV). In this case one 
can choose between the two opposite charged I¥-bosons, namely W~ (IWCH = 
1) and W + (IWCH = 2), and the invariant-mass of interest must be entered by 
the user in GeV. 

OWIDTH (CHARACTER* 1) 

It allows a different choice of the value of the VF-width. OWIDTH = Y means 
that the tree-level Standard Model formula for the VF-width is used; OWIDTH = 
N requires that the W- width is supplied by the user in GeV. This freedom in the 
choice of the VF-width has been introduced so as to simplify the comparisons 
of the results of our code with those of other programs. 

NSCH( INTEGER) 

The value of NSCH allows the user to choose the calculational scheme for the 
weak mixing angle and the gauge coupling. Three choices are available. If 
NSCH=1, the input parameters used are Gp, Myy, Mz and the calculation is 
performed at tree level. If NSCH = 2 or 3, the input parameters used are 
a(Q 2 ), Gf, Mw or a(Q 2 ), Gf, Mz, respectively, and the calculation is per- 
formed using the QED coupling constant at a proper scale Q 2 , which is re- 
quested as further input. 

OCOUL (CHARACTER* 1) 

This flag allows the user to switch on/off the contribution of the Coulomb 
correction. If OCOUL = Y (N) the distributions are computed with (without) 
the above electromagnetic effect. 

5 Test Run Output 

The typical calculations that can be performed with WWGENPV are illustrated 
by the following examples of test run output. 
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Sample 1 Calculation of total cross section, average energy loss and average 
invariant-mass loss for 10 8 weighted events, including initial-state QED 
corrections. 

Sample 2 Calculation of total cross section for 10 4 unweighted events and 
filling the n-tuple EXAMPLE . DAT, including initial-state QED corrections. 

The output of sample 1 can be obtained with the following input: 
OGEN = 'I' 

The integration branch is chosen. 
RS = 176. DO 

The centre-of-mass energy is fixed at 176 GeV. 

NHITWMAX = 100000000 

The maximum number of calls is fixed. 

I QED = 1 

Initial-state QED corrections are included. 
0DIS = 'T' 

The total cross section is computed (together with the corresponding average 
energy and invariant-mass loss). 

0WIDTH = 'Y' 

The W^-boson width is computed at the tree level in the Standard Model. 
NSCH = 2 

a(Q 2 ), Gf and Mw are used as input parameters. 

1/ALPHAI = 128.07D0 

The value of l/a(Q 2 ) is entered. 

0C0UL = 'N' 
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The Coulomb correction is not included. 

The output of sample 2 can be obtained with the following modifications: 
OGEN = 'G' 

The generation branch is chosen. 
NHITMAX = 10000 

The number of generated events is fixed. 
Moreover the value of ODIS is not required. 

Both for the integration and generation branch, the output returns the 
values of the "Feynman rules" GVE, GAE, GWF and GWWZ linked to the vertices 
occurring in the diagrams of e + e~ — > W + W~ — > 4/ as follows: 

Z f f vertex = 7^ (GVE + GAE 75) 
W u a M vertex = GWF 7m (1 + 75) 
WWZ vertex = GWWZ x Lorentz structure. 



6 Program Implementation 

In this section a short description of each routine implemented in the program 
is given. For random-number generation the routine RANLUX by F. James is 
used and inserted in the code. For more details about RANLUX the reader is 
referred to the proper literature 0. 

SUBROUTINE CUTUSER (EVENT , CUTFLAG) 
REAL* 8 CUTFLAG 
DIMENSION EVENT (15) 

To be supplied by the user. The subroutine allows the user to implement 
any kind of experimental cut on the four-momenta of the final-state fermions. 
Globally, the event is specified by 15 ingredients: the radiative variables x\ t i 
(EVENT (1) , EVENT (2)), the beam energy (EVENT (3)) and the components of 
the three-momenta of the four outgoing fermions (EVENT (I), with I — 4, 15 ). 
The energies are simply derived by the mass-shell condition for ultrarelativistic 
fermions. The reference frame is assumed to be with the z-axis along the 
direction of the incoming positron. 

The routine returns a flag (CUTFLAG = 0, 1), which controls if the event is 
accepted (CUTFLAG = 1) or rejected (CUTFLAG = 0) by the selection criteria. 
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SUBROUTINE BORN (S , X13 , X14 , X34 , X35 , X15 , X16 , X23 , X24 , 

X25 , X26 , X36 , X45 , X46 , X56 , SQM , IOVFL , IOEFL , 
IOES, JOEFL, JOES) 

IMPLICIT REAL* 8 (A-H.O-Z) 

It returns in SQM the squared tree-level matrix element associated to the s- and 
t-channels Feynman diagrams of the kernel process e + e~ — ► W + W~ Af. 
The squared matrix element is computed numerically, squaring the ampli- 
tude obtained by means of SCHOONSCHIP M in the framework of the helicity- 
amplitude formalism || for massless final fermions. If the effect of the Coulomb 
correction is required, the value returned for SQM is given by: 

SQM = SQM ■ {1 + a/v5c} , (16) 

where 5c represents the Coulomb correction factor and v the relative velocity 
of the two H^-bosons (see || for more details and explicit formulae). 

FUNCTION DOP(X) 

IMPLICIT REAL* 8 (A-H.O-Z) 

COMMON/EM/RL , ETA , BETA , SDELTA 



The QED structure function, adopted to account for the radiation emitted 
by the incoming leptons, is computed according to the analytical expression 
reported above. It is given by the sum of an exponentiated term of the Gribov- 
Lipatov form, associated to soft multiphoton emission, with up to 0(a 2 ) contri- 
butions describing hard bremsstrahlung in collinear approximation. Actually, 
in order to perform importance sampling on the generation of the radiative 
variable x in the soft region x — > 1, the following formula is implemented in 
FUNCTION DOP: 



D(x) = D(x)/p(x), (17) 



where p(x) is the weight function given by fll4D . 

Acknowledgements 

The authors would like to thank the INFN, Sezione di Pavia, for having 
provided computer resources. They are also grateful to Giorgio Fumagalli and 
Valerio Vercesi for helpful assistance in numerical and graphical problems. Use- 
ful discussions and exchanges of information with D. Bardin and D. Charlton 
are gratefully acknowledged. 



12 



References 

[1] M. Cacciari, A. Deandrea, G. Montagna and O. Nicrosini, Z. Phys. C52 
(1991) 421; 

D. Bardin, M. Bilenky, A. Olchevsky and T. Riemann, Phys. Lett. B308 
(1993) 403; 

V. Fadin, V.A. Khoze and A.D. Martin, Phys. Lett. B311 (1993) 311. 

[2] F. A. Berends, R. Kleiss and R. Pittau, Nucl. Phys. B424 (1994) 308 and 
B426 (1994) 344; Comput. Phys. Commun. 85 (1995) 437; 
H. Anlauf et al, Comput. Phys. Commun. 79 (1994) 487; 
G. J. van Oldenborgh, P. J. Franzini and A. Borrelli, Comput. Phys. 
Commun. 83 (1994) 14; 

T. SjSstrand, Comput. Phys. Commun. 82 (1994) 74; CERN-TH. 7112/93. 

[3] G. Montagna, O. Nicrosini, G. Passarino and F. Piccinini, Phys. Lett. 
B348 (1995) 178. 

[4] E. A. Kuraev and V. S. Fadin, Sov. J. Nucl. Phys. 41 (1985) 466; 

G. Altarelli and G. Martinelli, in Physics at LEP, CERN Report 86-02, 
J. Ellis and R. Peccei eds. (CERN, Geneva, 1986); 

O. Nicrosini and L. Trentadue, Phys. Lett. B196 (1987) 551; Z. Phys. C39 
(1988) 479. 

[5] F. James, Rep. Prog. Phys. 34 (1980) 1145. 

[6] V. S. Fadin, V. A. Khoze, A. D. Martin and A. Chapovsky, Coulomb Ef- 
fects in W + W~ Production, Preprint DTP/94/116 and references therein. 

[7] F. James, Comput. Phys. Commun. 79 (1994) 111. 

[8] SCH00NSCHIP, a program for symbol handling by M. Veltman, see 

H. Strubbe, Comput. Phys. Commun. 8 (1974) 1. 

[9] G. Passarino, Nucl. Phys. B237 (1984) 249. 



13 



TEST RUN OUTPUT 



Sample 1 



SQRT(S) = 176.0 GEV 

M_Z = 91.1888 GEV 

G_Z = 2.4974 GEV 

M_W = 80.23 GEV 

G_W = 2.03367033063195 GEV 

STH2 = 0.2310309124510679 

GVE = -1.409737267299689E-02 
GAE = -0.185794027211796 
GWF = 0.230409927395451 
GWWZ = -0.571479454308384 

RANLUX INITIALIZED BY RLUXGO FROM SEEDS 2 
NCALLS = 100000000 
XSECT FOR WEIGHTED EVENTS 
EFF = 0.999998000004 

XSECT = 13.52642598410257 +- 2 . 200641628422248E-03 (PB) 
ENERGY LOSS 

MEAN EG = 1.18064898022852 +- 6 . 041171111535717E-04 (GEV) 
INVARIANT-MASS LOSS 

MEAN IML = 1.17777558571609 +- 6 . 023399436481718E-04 (GEV) 



Sample 2 



SQRT(S) = 176.0 GEV 

M_Z = 91.1888 GEV 

G_Z = 2.4974 GEV 

M_W = 80.23 GEV 

G_W = 2.03367033063195 GEV 

STH2 = 0.2310309124510679 

GVE = -1.409737267299689E-02 
GAE = -0.185794027211796 
GWF = 0.230409927395451 
GWWZ = -0.571479454308384 

RANLUX INITIALIZED BY RLUXGO FROM SEEDS 2 
NHIT = 10000 

XSECT FOR UNWEIGHTED EVENTS 
EFF = 1.341451799624662E-02 

XSECT = 13.4145179962466 +- 0.133242395670936 (PB) 
NBIAS/NMAX = 0.0 



